/*************************************************************/
/*   Headers for a neo-hookean elasticity tensor in MOOSE    */
/*                  Status: not working                      */
/*                                                           */
/* Copyright (C) 2018  eDgar                                 */
/*                                                           */
/* This program is free software: you can redistribute it    */
/* and/or modify it under the terms of the GNU General       */
/* Public License as published by the Free Software          */
/* Foundation, either version 3 of the License, or (at your  */
/* option) any later version.                                */
/*                                                           */
/* This program is distributed in the hope that it will be   */
/* useful, but WITHOUT ANY WARRANTY; without even the        */
/* implied warranty of MERCHANTABILITY or FITNESS FOR A      */
/* PARTICULAR PURPOSE.  See the GNU General Public License   */
/* for more details.                                         */
/*                                                           */
/* You should have received a copy of the GNU General Public */
/* License along with this program.  If not, see             */
/* <http://www.gnu.org/licenses/>.                           */
/*************************************************************/

// We need to get the specification of how to call the
// functions which will be defined below. That is the main
// point of the header (.h) file
#include "ComputeElasticityTensorNeohookean.h"

// Populate the template which was sketched in the header
// file. The parameters given here can be used by the user
// in the input file.
template <>
InputParameters
validParams<ComputeElasticityTensorNeohookean>()
{
  InputParameters params = validParams<ComputeElasticityTensorBase>();

  // Description which is shown in the auto-generated
  // documentation
  params.addClassDescription
    ("DO NOT USE!!\n"
     "ComputeElasticityTensorNeohookean defines a "
     "elasticity tensor for the hyperelastic"
     "neo-Hookean material. "
     "\\["
     "c_{ijkl} = \\frac{\\lambda_{0}}{J}"
     " \\cdot\\delta_{ij}\\cdot\\delta_{kl}"
     " + \\frac{"
     "    \\mu_{0} - \\lambda_{0} \\cdot \\ln\\left("
     "      J\\right)}{J} \\cdot \\left("
     "        \\delta_{ik}\\cdot\\delta_{jl}"
     "        + \\delta_{il}\\cdot\\delta_{kj}"
     "      \\right)"
     "\\]"
     );

  // Give a default value to \mu (set -1; if the value is
  // not set to a positive value by the user, there is a
  // check below which throws an error)
  params.addParam<Real>
    ("mu",
     -1.0,
     "Lamé second constant (shear) of the material.");
  // Allow lambda to be an input parameter
  params.addParam<Real>
    ("lambda",
     "Lame's first constant (non-shear) of the material.");

  // Allow Poisson's ratio and Young's modulus as input
  // parameters
  params.addParam<Real>
    ("nu", "Poisson's ratio of the material.");
  params.addParam<Real>
    ("E", "Young's modulus of the material.");

  return params;
}

// Notice that the parameters become constants
// when we make this constructor
ComputeElasticityTensorNeohookean::
ComputeElasticityTensorNeohookean
(const InputParameters & parameters)
  : ComputeElasticityTensorBase(parameters),

    // When this function is called, it is because the
    // application is compiled and executed with an input
    // file (.i).
    //
    // First, make sure that the parameters were given in
    // the .i (input) file. Note that we called it "lambda"
    // in the template above, but the variable in the .h
    // file is _lambda. For the user filling the .i (input)
    // file, it will be `lambda'. For calculations in the
    // code, it will be `_lambda'.
    _lambda_set(parameters.isParamValid("lambda")),
    _mu_set(parameters.isParamValid("mu")),
    _nu_set(parameters.isParamValid("nu")),
    _E_set(parameters.isParamValid("E")),

    // Get parameters from the input file. The following
    // line is an implicit conditional. It checks if the
    // Boolean value is True, in which case, it assigns the
    // value after the `?' (and before `:'). If the value is
    // False, then the variable gets the value indicated
    // after `:'. In this case, if `lambda' was set in the
    // .i (input) file, it gets the value from
    // there. Otherwise, it sets it to -1.
    _lambda(_lambda_set ? getParam<Real>("lambda") : -1.0),
    // The same for `mu'
    _mu(_mu_set ? getParam<Real>("mu") : -1.0),

    // `nu' and `E' are parameters, because I don't need to
    // store them for anything else (so far).
    _nu(_nu_set ? getParam<Real>("nu") : -1.0),
    _E(_E_set ? getParam<Real>("E") : -1.0),

    // Take the deformation gradient from
    // ComputeFiniteStrain
    _deformation_gradient(getMaterialProperty<RankTwoTensor>
                          ("deformation_gradient")),

    // "You can declare material property for J_1 in the
    // ComputeElasticityTensorNeohookean.C and get material
    // property J_1 in ComputeStressNeohookean.C." --Wen Jiang
    //
    // I hope that this is the way to declare a new property
    // (based on FiniteStrainCrystalPlasticity)
    _J_1(declareProperty<Real>("jacobian_inverse")),
    _logJ(declareProperty<Real>("jacobian_log"))
{
  // Neo-hookean material model

  // Now we are going to initialise variables for the analysis
  // TODO:

  // Make sure that only one of these are given:
  // 1. Young's modulus and the Poisson's ratio
  // 2. Lamé's constants lambda and mu
  if (! ( (_E_set && _nu_set) || (_lambda_set && _mu_set) ) )
    mooseError
      ("One of these combinations is needed:\n"
       "(1) Young's modulus (E) and the Poisson's ratio (nu)\n"
       "(2) Lamé's constants lambda and mu\n"
       "in the input file for the neo-Hoookean\n"
       "'" + name() + "'.");
  else if (_E_set && _nu_set) {
      // If both Young's modulus and Poisson's ratio are
      // given, use them to calculate Lamé's constants
      if (_nu_set && (_nu <= -1.0 || _nu >= 0.5))
        mooseError
          ("Poissons ratio must be greater than -1.0 and less"
           " than 0.5 in material '" + name() + "'.");
      if (_E_set && _E <= 0.0)
        mooseError
          ("Youngs modulus must be positive in material '"
           + name() + "'.");

      // Value of Lamé constants in terms of Young's modulus
      // and Poisson's ratio
      _mu = _E / (2.0 * (_nu + 1.0));
      // // 3 multiplications, 1 division, 1 addition, 1 subtraction:
      // _lambda = E * _nu/ ((_nu + 1) * (1 - 2 * _nu));
      // 3 multiplications, 1 division, 1 subtraction, less clear:
      _lambda = 2.0 * _mu * _nu / (1.0 - 2.0 * _nu);
    }
  else if (_lambda_set && _mu_set) {
    // Make sure that the values of Lamé's constants are
    // not both negative or zero.
    // Holzapfel, G. A., Nonlinear solid mechanics : a
    // continuum approach for engineering (2000), :
    // Wiley. p.250
    if ((_lambda < 0.0 || _mu < 0.0) ||
        (_lambda == 0.0 && _mu == 0.0))
      mooseError("The Lamé constants must be positive in material '"
                 + name() + "'.");
  }

  // All tensors created by this class are always isotropic
  // (ComputeFiniteStrainElasticStress can only be used with
  // elasticity tensors that guarantee isotropy).
  // TODO: how do I fill it in _and_ check if it's isotropic?
  /* if (_c_ijkl.isSymmetric()) */
  issueGuarantee(_elasticity_tensor_name, Guarantee::ISOTROPIC);
}

void
ComputeElasticityTensorNeohookean::computeQpElasticityTensor()
{
  // The elasticity tensor of the neo-Hookean model
  // in Eulerian coordinates is
  //
  // LaTeX:
  // c_{ijkl} = \frac{\lambda}{J}
  //                \cdot\delta_{ij}\cdot\delta_{kl}
  //                + \frac{
  //                   \mu - \lambda \cdot \ln\left(
  //                     J\right)}{J} \cdot \left(
  //                       \delta_{ik}\cdot\delta_{jl}
  //                       + \delta_{il}\cdot\delta_{kj}
  //               \right)
  //
  // - Belytschko, T., Liu, W. K., Moran, B., & Elkhodary,
  //   K. (2014). Nonlinear finite elements for continua and
  //   structures. : John Wiley and Sons, Ltd.

  // You have to make sure that whatever operatations
  // (including addition, subtraction, etc. work as you
  // think). In this case, check that the definition of
  // RankFourTensor has the required operations.
  // [[moose/framework/src/utils/RankFourTensor.C]]

  // LaTeX:
  // \( \left(\delta_{ik}\,\delta_{jl}
  //    + \delta_{il}\,\delta_{jk}\right) \)
  //
  // 1.0  0.0  0.0  0.0  0.0  0.0
  // 0.0  1.0  0.0  0.0  0.0  0.0
  // 0.0  0.0  1.0  0.0  0.0  0.0
  // 0.0  0.0  0.0  0.5  0.0  0.0
  // 0.0  0.0  0.0  0.0  0.5  0.0
  // 0.0  0.0  0.0  0.0  0.0  0.5
  // initIdentitySymmetricFour() in
  // [[moose/framework/src/utils/RankFourTensor.C]]
  //
  // initIdentitySymmetricFour
  // gives (de_ij de_kl + de_il de_jk)/2 -- andrew.wilkins
  //
  // It's possible that someone might use this with crazy
  // non-symmetric stresses and strains and
  // initIdentitySymmetricFour essentially symmetrises the
  // strain before forming the stress. -- andrew.wilkins
  //
  // See also [BelytschkoLiuMoranElkhodary14, p.240]
  const RankFourTensor dik_djl__dil_dkj = RankFourTensor::initIdentitySymmetricFour;

  // Identity matrix
  // 1  0  0
  // 0  1  0
  // 0  0  1
  // initIdentity() in
  // [[moose/framework/src/utils/RankTwoTensor.C]]
  const RankTwoTensor dij(RankTwoTensor::initIdentity);

  // LaTeX:
  // \( \delta_{ij} \, \delta_{kl} \)
  //
  // 1  1  1  0  0  0
  // 1  1  1  0  0  0
  // 1  1  1  0  0  0
  // 0  0  0  0  0  0
  // 0  0  0  0  0  0
  // 0  0  0  0  0  0
  const RankFourTensor dij_dkl = dij.outerProduct(dij);

  // // Another way for dij_dkl
  // // To make sure that nothing is out of place, fill the
  // // whole tensor with zeros
  // // [[moose/framework/include/utils/RankFourTensor.h]]
  // dij_dkl.zero();
  // // Set the corresponding values
  // // [[moose/framework/src/utils/RankFourTensor.C]]
  // dij_dkl.fillGeneralIsotropic(1, 0, 0);

  // Deformation gradient tensor F = 1 + du_{i}/dX_{j}
  RankTwoTensor F;
  // "You can declare material property for J_1 in the
  // ComputeElasticityTensorNeohookean.C and get material
  // property J_1 in ComputeStressNeohookean.C." --Wen Jiang
  // Determinant of the deformation gradient tensor
  // \(\det{F_{ij}}\)
  //
  // Instead of this:
  // Real _J = _deformation_gradient[_qp].det();
  // Do this:
  // Jacobian (determinant of the deformation gradient)
  Real J;

  // Just copy the value for the given quadrature point
  F = _deformation_gradient[_qp];
  // Calculate Jacobian with the determinant
  J = F.det();
  // Calculate inverse of Jacobian
  _J_1[_qp] = 1.0 / J;
  // Calculate log(J)
  _logJ[_qp] = log(J);

  // Knowing that the * operator is defined as we expect
  // between a tensor and a scalar, as well as the +
  // operator between the two RankFourTensor
  // [[moose/framework/src/utils/RankFourTensor.h]]
  // , from the way the neo-hookean material is
  // defined:
  _c_ijkl = _J_1[_qp] * (_lambda * dij_dkl
                   + (_mu - _lambda * _logJ[_qp]) * 2.0 * dik_djl__dil_dkj);

  // The latter is equivalent to this. Seems different to
  // Bonet & Wood eq.5.37, p.149/279:
  // // Push forward of the fourth order elasticity tensor
  // // Green-Cauchy right deformation tensor
  // RankTwoTensor Cij = F.transpose() * F;
  // // Inverse of Cij
  // RankTwoTensor Cij_1 = Cij.inverse();
  // // Elasticity tensor in Lagrange coordinates
  // RankFourTensor C_ijkl;
  // // Fourth order identity tensor I_ijkl = Cik^(-1.0) * Cjl^(-1.0)
  // RankFourTensor I_ijkl;
  //
  // // Fdik_djl: Push-forward of I_ijkl
  // RankFourTensor Fdik_djl;
  // // \LaTeX:
  // // B_{ijkl} = F_{iI} F_{jJ} F_{kK} F_{lL} A_{IJKL}
  //
  // I_ijkl.zero();
  // for (int i = 0; i < 3; ++i)
  //   {
  //     for (int j = 0; j < 3; ++j)
  //       {
  //         for (int k = 0; k < 3; ++k)
  //           {
  //             for (int l = 0; l < 3; ++l)
  //               {
  //                 I_ijkl(i, j, k, l) += Cij_1(i, k) * Cij_1(j, l);
  //               }
  //           }
  //       }
  //   }
  //
  // // "There is a rotate method in RankFourTensor.C that
  // // performs the rank-2 and rank-4 operation you've listed."
  // // -- Stephanie Pitts
  // Fdik_djl = I_ijkl;
  // Fdik_djl.rotate(F);
  // _c_ijkl = _J_1[_qp] * (_lambda * dij_dkl
  //                  + 2 * (_mu - _lambda * _logJ) * Fdik_djl);
  //

  // This assigns an elasticity tensor at each quadrature
  // point (_qp) to each of the elements of the material.
  _elasticity_tensor[_qp] = _c_ijkl;
}
